Exp I - microbiome community composition
19 August 2024
1 Setup
1.1 Libs
## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-6.1library("phyloseq")
library("ggplot2")
library("cowplot")
# #install.packages("colorBlindness")
# library("colorBlindness")
#BiocManager::install("microbiome")
#library("microbiome")
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union#install.packages("microViz", repos = c(davidbarnett = "https://david-barnett.r-universe.dev", getOption("repos")))
library(microViz)## microViz version 0.12.4 - Copyright (C) 2021-2024 David Barnett
## ! Website: https://david-barnett.github.io/microViz
## ✔ Useful? For citation details, run: `citation("microViz")`
## ✖ Silence? `suppressPackageStartupMessages(library(microViz))`##
## Attaching package: 'ggpubr'## The following object is masked from 'package:microViz':
##
## stat_chull## The following object is masked from 'package:cowplot':
##
## get_legend##
## Attaching package: 'speedyseq'## The following objects are masked from 'package:phyloseq':
##
## filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
## tip_glom, transform_sample_counts1.2 Data
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 161 taxa and 458 samples ]:
## sample_data() Sample Data: [ 458 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columnsps.mq1 <- subset_samples(ps.clean,type=="A.albopictus")
ps.mq <- prune_taxa(taxa_sums(ps.mq1) > 0, ps.mq1)
ps.mq## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 89 taxa and 195 samples ]:
## sample_data() Sample Data: [ 195 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 89 taxa by 11 taxonomic ranks ]:
## taxa are columnsps.mw1 <- subset_samples(ps.clean,type=="Microbial Water")
ps.mw <- prune_taxa(taxa_sums(ps.mw1) > 0, ps.mw1)
ps.mw## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 49 taxa and 211 samples ]:
## sample_data() Sample Data: [ 211 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 49 taxa by 11 taxonomic ranks ]:
## taxa are columns##datasets but trimmed for mixed model things
ps.trim.mq <- readRDS("../04c.mixed_models/ps.trim.mq.rds")
ps.trim.mq #29 taxa, 195 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 29 taxa and 195 samples ]:
## sample_data() Sample Data: [ 195 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 29 taxa by 11 taxonomic ranks ]:
## taxa are columns## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 39 taxa and 211 samples ]:
## sample_data() Sample Data: [ 211 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 39 taxa by 11 taxonomic ranks ]:
## taxa are columns2 PCOA
2.1 All samples
By ‘all’ here I really mean without separating the two experimental types (mosquitoes & mesocosm water)
2.1.1 Relative abundance & Bray
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 161 taxa and 458 samples ]:
## sample_data() Sample Data: [ 458 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columnsps.all.rel1 <- subset_samples(ps.clean.rel,type=="A.albopictus"|type=="Microbial Water")
ps.all.rel1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 161 taxa and 406 samples ]:
## sample_data() Sample Data: [ 406 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columns## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 103 taxa and 406 samples ]:
## sample_data() Sample Data: [ 406 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 103 taxa by 11 taxonomic ranks ]:
## taxa are columnsord.all.rel <- ordinate(ps.all.rel, "PCoA", "bray")
plot.all.rel <- plot_ordination(ps.all.rel, ord.all.rel,color="type")+
stat_ellipse()+
theme_cowplot()
plot.all.rel #beautiful U shape
2.1.2 Aitchison (spelling?)
2.1.2.1 Full dataset
ps.clean.exp <- subset_samples(ps.clean,type=="A.albopictus"|type=="Microbial Water")
ps.exp.clr <- tax_transform(ps.clean.exp,"clr")
ord.exp.clr <- ordinate(ps.exp.clr, "PCoA", "euclidean")
gg.exp <- plot_ordination(ps.exp.clr, ord.exp.clr,color="type",label="mesocosm",shape="temperature")+
scale_color_manual(name="Type",values=c("grey20","grey50"),labels=c("Mosq.","Water"))+
theme_cowplot()+
scale_shape_manual(name="Temperature",values=c(15,17),labels=c("Cool","Warm"))+
ggtitle("Full dataset")
#scale_fill_manual(values=c("white","black"),name="Type",labels=c("Mosq.","Water"))
gg.exp
##extract data so I can make it look how I want
df.exp.plot <- gg.exp[["data"]]
##new variable for plotting colors
df.exp.plot$inf.temp <- paste0(df.exp.plot$infusion,"_",df.exp.plot$temperature)
gg.exp.full <- ggplot(df.exp.plot,aes(x=Axis.1,y=Axis.2,fill=inf.temp,shape=inf.temp))+
theme_cowplot()+
stat_ellipse(aes(linetype=type,color=inf.temp),linewidth=0.6)+
xlab("Axis 1 (52.3%)")+
ylab("Axis 2 (11.4%)")+
geom_point(color="gray70")+
scale_color_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
ggtitle("Both sample types")+
scale_linetype_manual(values=c("solid","dashed"),name="Sample type",labels=c("Mosq.","Water"))+
theme(legend.key.width=unit(1,"cm"),plot.title = element_text(face = "plain"))
#theme(axis.text.x=element_text(angle=90),legend.key.width=unit(1,"cm"))
gg.exp.full
2.2 Mosquitoes - less trimmed
2.2.1 Multivariate plot
#ps.mq.rel <- transform_sample_counts(ps.trim.mq, function(x) x / sum(x) )
#ps.mq.rel <- transform_sample_counts(ps.mq, function(x) x / sum(x) )
#ps.trim.mq.hel <- transform_sample_counts(ps.trim.mq, function(OTU) sqrt(OTU/sum(OTU)))
ps.mq.clr <- tax_transform(ps.mq,"clr")
ord.mq.clr <- ordinate(ps.mq.clr, "PCoA", "euclidean")
gg.mq <- plot_ordination(ps.mq.clr, ord.mq.clr,color="infusion")+
stat_ellipse()+
theme_cowplot()
gg.mq
##extract data so I can make it look how I want
df.plot.mq <- gg.mq[["data"]]
##new variable for plotting colors
df.plot.mq$inf.temp <- paste0(df.plot.mq$infusion,"_",df.plot.mq$temperature)
gg.pcoa.aitch.mq <- ggplot(df.plot.mq,aes(x=Axis.1,y=Axis.2,fill=inf.temp,shape=inf.temp))+
theme_cowplot()+
stat_ellipse(aes(color=inf.temp),linewidth=0.6)+
xlab("Axis 1 (17.6%)")+
ylab("Axis 2 (14.3%)")+
geom_point(color="gray70")+
scale_color_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
ggtitle("Mosquitoes only")+
#scale_linetype_manual(values=c("solid","dashed","solid","dashed","solid","dashed"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
theme(legend.key.width=unit(1,"cm"),plot.title = element_text(face = "plain"))
#theme(axis.text.x=element_text(angle=90),legend.key.width=unit(1,"cm"))
gg.pcoa.aitch.mq
2.2.2 Both panels - with diversity (Fig. 2)
Not knitting because requires the diversity script to be run also
2.3 Comparing prevalence-trimmed levels
ps.trim.exp <- merge_phyloseq(ps.trim.mq,ps.trim.mw)
ps.trim.exp.clr <- tax_transform(ps.trim.exp,"clr")
ord.trim.exp.clr <- ordinate(ps.trim.exp.clr, "PCoA", "euclidean")
gg.exp.trim <- plot_ordination(ps.trim.exp.clr, ord.trim.exp.clr,color="type",label="mesocosm",shape="temperature")+
scale_color_manual(name="Type",values=c("grey20","grey50"),labels=c("Mosq.","Water"))+
theme_cowplot()+
scale_shape_manual(name="Temperature",values=c(15,17),labels=c("Cool","Warm"))+
ggtitle("Prevalence-trimmed")
#scale_fill_manual(values=c("white","black"),name="Type",labels=c("Mosq.","Water"))
gg.exp.trim
3 Trimmed vs. untrimmed stats
3.1 Mesocosm summed
3.1.1 Data - trimmed for mixed models
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionsamdf.trim.mq <- data.frame(sample_data(ps.trim.mq))
samdf.trim.mq.less <- samdf.trim.mq %>%
dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.trim.mq.uniq <- distinct(samdf.trim.mq.less)
row.names(samdf.trim.mq.uniq) <- samdf.trim.mq.uniq$mesocosm
sample_data(ps.trim.mq.mer) <- sample_data(samdf.trim.mq.uniq)
##merge by mesocosm - water
ps.trim.mw.mer <- merge_samples(ps.trim.mw, "mesocosm")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionsamdf.trim.mw <- data.frame(sample_data(ps.trim.mw))
samdf.trim.mw.less <- samdf.trim.mw %>%
dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.trim.mw.uniq <- distinct(samdf.trim.mw.less)
row.names(samdf.trim.mw.uniq) <- samdf.trim.mw.uniq$mesocosm
sample_data(ps.trim.mw.mer) <- sample_data(samdf.trim.mw.uniq)3.1.2 Statz trimmed - mosquitoes
ps.trim.mq.mer.clr <- tax_transform(ps.trim.mq.mer,"clr")
seq.trim.mq.mer.clr <- data.frame(ps.trim.mq.mer.clr@otu_table)
samdf.trim.mq.mer.clr <- data.frame(ps.trim.mq.mer.clr@sam_data)
dist.aich.trim.mq.mer <- vegdist(seq.trim.mq.mer.clr, method="euclidean")
##beta - infusion
bet.aich.trim.mq.mer.inf <- betadisper(dist.aich.trim.mq.mer,samdf.trim.mq.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mq.mer.inf,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 10.345 5.1723 0.9967 999 0.372
## Residuals 55 285.407 5.1892##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 10.345 5.1723 0.9967 999 0.377
## Residuals 55 285.407 5.1892
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## OL SG SW
## OL 0.92400 0.151
## SG 0.92100 0.314
## SW 0.16761 0.33496
##beta - temperature
bet.aich.trim.mq.mer.tem <- betadisper(dist.aich.trim.mq.mer,samdf.trim.mq.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mq.mer.tem,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 3.026 3.0261 0.5443 999 0.462
## Residuals 56 311.330 5.5595
##beta - dispersal
bet.aich.trim.mq.mer.dis <- betadisper(dist.aich.trim.mq.mer,samdf.trim.mq.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mq.mer.dis,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 4.359 4.3591 0.8315 999 0.35
## Residuals 56 293.560 5.2421
##adonis
adonis2(dist.aich.trim.mq.mer ~ infusion+temperature+dispersal, data=samdf.trim.mq.mer.clr, permutations=999) ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.aich.trim.mq.mer ~ infusion + temperature + dispersal, data = samdf.trim.mq.mer.clr, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## infusion 2 759.9 0.15531 5.1886 0.001 ***
## temperature 1 133.2 0.02723 1.8193 0.045 *
## dispersal 1 118.7 0.02426 1.6213 0.063 .
## Residual 53 3881.3 0.79320
## Total 57 4893.2 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.3 Statz trimmed - water
ps.trim.mw.mer.clr <- tax_transform(ps.trim.mw.mer,"clr")
seq.trim.mw.mer.clr <- data.frame(ps.trim.mw.mer.clr@otu_table)
samdf.trim.mw.mer.clr <- data.frame(ps.trim.mw.mer.clr@sam_data)
dist.aich.trim.mw.mer <- vegdist(seq.trim.mw.mer.clr, method="euclidean")
##beta - infusion
bet.aich.trim.mw.mer.inf <- betadisper(dist.aich.trim.mw.mer,samdf.trim.mw.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mw.mer.inf,permutations=999) #sig***##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 73.567 36.784 22.286 999 0.001 ***
## Residuals 69 113.885 1.651
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 73.567 36.784 22.286 999 0.001 ***
## Residuals 69 113.885 1.651
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## OL SG SW
## OL 1.0000e-03 0.013
## SG 4.1802e-04 0.001
## SW 1.1631e-02 3.8675e-08
##beta - temperature
bet.aich.trim.mw.mer.tem <- betadisper(dist.aich.trim.mw.mer,samdf.trim.mw.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mw.mer.tem,permutations=999) #0.057 .##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 12.387 12.3865 3.8437 999 0.051 .
## Residuals 70 225.580 3.2226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##beta - dispersal
bet.aich.trim.mw.mer.dis <- betadisper(dist.aich.trim.mw.mer,samdf.trim.mw.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.trim.mw.mer.dis,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 8.585 8.5847 2.5233 999 0.127
## Residuals 70 238.149 3.4021
##adonis
adonis2(dist.aich.trim.mw.mer ~ infusion+temperature+dispersal, data=samdf.trim.mw.mer.clr, permutations=999) ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.aich.trim.mw.mer ~ infusion + temperature + dispersal, data = samdf.trim.mw.mer.clr, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## infusion 2 5766.5 0.64480 67.6368 0.001 ***
## temperature 1 213.1 0.02382 4.9982 0.004 **
## dispersal 1 107.4 0.01201 2.5196 0.034 *
## Residual 67 2856.1 0.31936
## Total 71 8943.1 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Df SumOfSqs R2 F Pr(>F)
# infusion 2 5766.5 0.64480 67.6368 0.001 ***
# temperature 1 213.1 0.02382 4.9982 0.002 **
# dispersal 1 107.4 0.01201 2.5196 0.053 .
# Residual 67 2856.1 0.31936
# Total 71 8943.1 1.00000
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 13.1.4 Data - not trimmed
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionsamdf.mq <- data.frame(sample_data(ps.mq))
samdf.mq.less <- samdf.mq %>%
dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.mq.uniq <- distinct(samdf.mq.less)
row.names(samdf.mq.uniq) <- samdf.mq.uniq$mesocosm
sample_data(ps.mq.mer) <- sample_data(samdf.mq.uniq)
##merge by mesocosm - water
ps.mw.mer <- merge_samples(ps.mw, "mesocosm")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion3.1.5 Statz not trimmed - mosquitoes
ps.mq.mer.clr <- tax_transform(ps.mq.mer,"clr")
seq.mq.mer.clr <- data.frame(ps.mq.mer.clr@otu_table)
samdf.mq.mer.clr <- data.frame(ps.mq.mer.clr@sam_data)
dist.aich.mq.mer <- vegdist(seq.mq.mer.clr, method="euclidean")
##beta - infusion
bet.aich.mq.mer.inf <- betadisper(dist.aich.mq.mer,samdf.mq.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.mq.mer.inf,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 10.42 5.2081 0.6266 999 0.571
## Residuals 55 457.13 8.3114##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 10.42 5.2081 0.6266 999 0.546
## Residuals 55 457.13 8.3114
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## OL SG SW
## OL 0.48000 0.574
## SG 0.48625 0.341
## SW 0.55470 0.32896
##beta - temperature
bet.aich.mq.mer.tem <- betadisper(dist.aich.mq.mer,samdf.mq.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.mq.mer.tem,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01 0.0080 0.001 999 0.978
## Residuals 56 438.96 7.8385
##beta - dispersal
bet.aich.mq.mer.dis <- betadisper(dist.aich.mq.mer,samdf.mq.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.mq.mer.dis,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.44 0.4433 0.0583 999 0.801
## Residuals 56 425.79 7.6034
##adonis
adonis2(dist.aich.mq.mer ~ infusion+temperature+dispersal, data=samdf.mq.mer.clr, permutations=999) ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.aich.mq.mer ~ infusion + temperature + dispersal, data = samdf.mq.mer.clr, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## infusion 2 888.5 0.11387 3.6025 0.001 ***
## temperature 1 206.8 0.02651 1.6772 0.027 *
## dispersal 1 171.6 0.02200 1.3918 0.113
## Residual 53 6535.6 0.83763
## Total 57 7802.5 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.6 Statz not trimmed - water
ps.mw.mer.clr <- tax_transform(ps.mw.mer,"clr")
seq.mw.mer.clr <- data.frame(ps.mw.mer.clr@otu_table)
samdf.mw.mer.clr <- data.frame(ps.mw.mer.clr@sam_data)
dist.aich.mw.mer <- vegdist(seq.mw.mer.clr, method="euclidean")
##beta - infusion
bet.aich.mw.mer.inf <- betadisper(dist.aich.mw.mer,samdf.mw.mer.clr$infusion)
#anova(bet.rad.uk2i)
permutest(bet.aich.mw.mer.inf,permutations=999) #sig***##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 81.949 40.974 32.543 999 0.001 ***
## Residuals 69 86.876 1.259
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 81.949 40.974 32.543 999 0.001 ***
## Residuals 69 86.876 1.259
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## OL SG SW
## OL 1.0000e-03 0.003
## SG 3.2943e-05 0.001
## SW 2.7862e-03 2.8516e-10
##beta - temperature
bet.aich.mw.mer.tem <- betadisper(dist.aich.mw.mer,samdf.mw.mer.clr$temperature)
#anova(bet.rad.uk2i)
permutest(bet.aich.mw.mer.tem,permutations=999) #0.057 .##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 13.948 13.9477 4.7498 999 0.03 *
## Residuals 70 205.555 2.9365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##beta - dispersal
bet.aich.mw.mer.dis <- betadisper(dist.aich.mw.mer,samdf.mw.mer.clr$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.aich.mw.mer.dis,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 7.431 7.4314 2.3515 999 0.142
## Residuals 70 221.221 3.1603
##adonis
adonis2(dist.aich.mw.mer ~ infusion+temperature+dispersal, data=samdf.mw.mer.clr, permutations=999) ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.aich.mw.mer ~ infusion + temperature + dispersal, data = samdf.mw.mer.clr, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## infusion 2 6047.0 0.62704 62.3525 0.001 ***
## temperature 1 233.9 0.02426 4.8240 0.007 **
## dispersal 1 113.9 0.01181 2.3491 0.051 .
## Residual 67 3248.9 0.33689
## Total 71 9643.7 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# adonis2(formula = dist.aich.mw.mer ~ infusion + temperature + dispersal, data = samdf.mw.mer.clr, permutations = 999)
# Df SumOfSqs R2 F Pr(>F)
# infusion 2 6353.4 0.63122 63.3973 0.001 ***
# temperature 1 240.6 0.02391 4.8023 0.005 **
# dispersal 1 114.0 0.01133 2.2749 0.053 .
# Residual 67 3357.2 0.33355
# Total 71 10065.2 1.00000
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 14 Unifrac things
4.1 Making the phylo tree info
Ran the following once & not doing it while knitting because it takes time
4.1.1 Subset needed taxa
Only need the taxa that end up in mosquitoes & mesocosm water samples
library("dada2")
library("DECIPHER")
library("phangorn")
setwd("~/Library/CloudStorage/GoogleDrive-nicolagk@hawaii.edu/My Drive/Mosquito_business/Mosquito_microbes/04.microbiome_analysis/04b.comm_comp")
ps.trim <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.trim #161 taxa, 458 samples
ps.exp <- subset_samples(ps.trim,type=="A.albopictus"|type=="Microbial Water")
ps.exp #406 samples
ps.exp.no0 <- prune_taxa(taxa_sums(ps.exp)!=0,ps.exp)
ps.exp.no0 #103 taxa4.1.2 Alignment & tree things
Tutorial from dada2 author here
tax.exp <- data.frame(tax_table(ps.exp.no0))
head(tax.exp)
seqs <- getSequences(tax.exp$Sequence)
names(seqs) <- row.names(tax.exp) # This propagates to the tip labels of the tree
names(seqs)
alignment <- AlignSeqs(DNAStringSet(seqs), anchor=NA)
phang.align <- phyDat(as(alignment, "matrix"), type="DNA")
dm <- dist.ml(phang.align)
treeNJ <- NJ(dm) # Note, tip order != sequence order
fit = pml(treeNJ, data=phang.align)
## negative edges length changed to 0!
fitGTR <- update(fit, k=4, inv=0.2)
fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
rearrangement = "stochastic", control = pml.control(trace = 0))
detach("package:phangorn", unload=TRUE)
#saveRDS(fitGTR, file="phylo.fit.all.rds")4.2 Unifrac objects
Can skip above stuff
setwd("~/Library/CloudStorage/GoogleDrive-nicolagk@hawaii.edu/My Drive/Mosquito_business/mosquito_microbes/04.microbiome_analysis/04b.comm_comp")
ps.tree <- readRDS("ps.exp.tree.rds")
ps.tree## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 103 taxa and 406 samples ]:
## sample_data() Sample Data: [ 406 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 103 taxa by 11 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 103 tips and 101 internal nodes ]:
## taxa are columns## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 29 taxa and 195 samples ]:
## sample_data() Sample Data: [ 195 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 29 taxa by 11 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 103 tips and 101 internal nodes ]:
## taxa are columns## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 39 taxa and 211 samples ]:
## sample_data() Sample Data: [ 211 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 39 taxa by 11 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 103 tips and 101 internal nodes ]:
## taxa are columns4.3 Mesocosm summed
4.3.1 Data - trimmed for mixed models
ps.trim.mq.rel <- transform_sample_counts(ps.trim.mq, function(x) x / sum(x))
ps.trim.mw.rel <- transform_sample_counts(ps.trim.mw, function(x) x / sum(x))
##merge by mesocosm - mosquitoes
ps.mq.mer <- merge_samples(ps.trim.mq.rel, "mesocosm")## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionsamdf.mq <- data.frame(sample_data(ps.trim.mq.rel))
samdf.mq.less <- samdf.mq %>%
dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.mq.uniq <- distinct(samdf.mq.less)
row.names(samdf.mq.uniq) <- samdf.mq.uniq$mesocosm
sample_data(ps.mq.mer) <- sample_data(samdf.mq.uniq)
ps.mq.mer## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 29 taxa and 58 samples ]:
## sample_data() Sample Data: [ 58 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 29 taxa by 11 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 29 tips and 28 internal nodes ]:
## taxa are columns## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercion
## Warning in asMethod(object): NAs introduced by coercionsamdf.mw <- data.frame(sample_data(ps.trim.mw.rel))
samdf.mw.less <- samdf.mw %>%
dplyr::select(mesocosm,infusion,dispersal,temperature)
samdf.mw.uniq <- distinct(samdf.mw.less)
row.names(samdf.mw.uniq) <- samdf.mw.uniq$mesocosm
sample_data(ps.mw.mer) <- sample_data(samdf.mw.uniq)
ps.mw.mer## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 39 taxa and 72 samples ]:
## sample_data() Sample Data: [ 72 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 39 taxa by 11 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 39 tips and 37 internal nodes ]:
## taxa are columns4.4 Statz mesocosm summed - mosquitoes
# Transform counts to relative abundances
ps.mq.mer.rel <- transform_sample_counts(ps.mq.mer, function(x) x / sum(x))
ps.mq.mer.rel## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 29 taxa and 58 samples ]:
## sample_data() Sample Data: [ 58 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 29 taxa by 11 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 29 tips and 28 internal nodes ]:
## taxa are columns#seq.mq.mer.rel <- data.frame(ps.mq.mer.rel@otu_table)
sam.mq.mer.rel <- data.frame(ps.mq.mer.rel@sam_data)
dist.wuni.mq.rel <- phyloseq::distance(ps.mq.mer.rel, method = "wunifrac")
##beta - infusion
bet.wuni.mq.rel.inf <- betadisper(dist.wuni.mq.rel,sam.mq.mer.rel$infusion)## Warning in betadisper(dist.wuni.mq.rel, sam.mq.mer.rel$infusion): some squared
## distances are negative and changed to zero##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.20598 0.102992 5.283 999 0.007 **
## Residuals 55 1.07223 0.019495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.20598 0.102992 5.283 999 0.006 **
## Residuals 55 1.07223 0.019495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## OL SG SW
## OL 0.7700000 0.002
## SG 0.7319054 0.028
## SW 0.0097278 0.0479558
##beta - temperature
bet.wuni.mq.rel.tem <- betadisper(dist.wuni.mq.rel,sam.mq.mer.rel$temperature)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mq.rel.tem,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.07522 0.075224 2.7138 999 0.101
## Residuals 56 1.55230 0.027720
## Warning in betadisper(dist.wuni.mq.rel, sam.mq.mer.rel$dispersal): some squared
## distances are negative and changed to zero##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.07251 0.072508 2.6287 999 0.106
## Residuals 56 1.54463 0.027583
##adonis
adonis2(dist.wuni.mq.rel ~ infusion+temperature+dispersal, data=sam.mq.mer.rel, permutations=999) ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.wuni.mq.rel ~ infusion + temperature + dispersal, data = sam.mq.mer.rel, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## infusion 2 0.42030 0.21125 7.7307 0.001 ***
## temperature 1 0.08700 0.04373 3.2005 0.049 *
## dispersal 1 0.04151 0.02086 1.5270 0.229
## Residual 53 1.44075 0.72415
## Total 57 1.98956 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Df SumOfSqs R2 F Pr(>F)
# infusion 2 0.42030 0.21125 7.7307 0.001 ***
# temperature 1 0.08700 0.04373 3.2005 0.039 *
# dispersal 1 0.04151 0.02086 1.5270 0.240
# Residual 53 1.44075 0.72415
# Total 57 1.98956 1.00000
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 14.5 Statz mesocosm summed - water
# Transform counts to relative abundances
ps.mw.mer.rel <- transform_sample_counts(ps.mw.mer, function(x) x / sum(x))
ps.mw.mer.rel## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 39 taxa and 72 samples ]:
## sample_data() Sample Data: [ 72 samples by 4 sample variables ]:
## tax_table() Taxonomy Table: [ 39 taxa by 11 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 39 tips and 37 internal nodes ]:
## taxa are columns#seq.mw.mer.rel <- data.frame(ps.mw.mer.rel@otu_table)
sam.mw.mer.rel <- data.frame(ps.mw.mer.rel@sam_data)
dist.wuni.mw.rel <- phyloseq::distance(ps.mw.mer.rel, method = "wunifrac")## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## Otu0049 -- in the phylogenetic tree in the data you provided.##beta - infusion
bet.wuni.mw.rel.inf <- betadisper(dist.wuni.mw.rel,sam.mw.mer.rel$infusion)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mw.rel.inf,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.000202 0.00010112 0.0716 999 0.933
## Residuals 69 0.097406 0.00141168##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.000202 0.00010112 0.0716 999 0.944
## Residuals 69 0.097406 0.00141168
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## OL SG SW
## OL 0.85900 0.695
## SG 0.83196 0.889
## SW 0.70404 0.87936
##beta - temperature
bet.wuni.mw.rel.tem <- betadisper(dist.wuni.mw.rel,sam.mw.mer.rel$temperature)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mw.rel.tem,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00373 0.0037270 0.7637 999 0.379
## Residuals 70 0.34161 0.0048801
##beta - dispersal
bet.wuni.mw.rel.dis <- betadisper(dist.wuni.mw.rel,sam.mw.mer.rel$dispersal)
#anova(bet.rad.uk2i)
permutest(bet.wuni.mw.rel.dis,permutations=999) #ns##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00071 0.0007117 0.1412 999 0.694
## Residuals 70 0.35294 0.0050421
##adonis
adonis2(dist.wuni.mw.rel ~ infusion+temperature+dispersal, data=sam.mw.mer.rel, permutations=999) ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.wuni.mw.rel ~ infusion + temperature + dispersal, data = sam.mw.mer.rel, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## infusion 2 1.04644 0.77043 168.1578 0.001 ***
## temperature 1 0.10069 0.07413 32.3606 0.001 ***
## dispersal 1 0.00265 0.00195 0.8519 0.409
## Residual 67 0.20847 0.15348
## Total 71 1.35825 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Df SumOfSqs R2 F Pr(>F)
# infusion 2 3.7520 0.78212 168.4687 0.001 ***
# temperature 1 0.2880 0.06004 25.8669 0.001 ***
# dispersal 1 0.0111 0.00231 0.9955 0.330
# Residual 67 0.7461 0.15552
# Total 71 4.7972 1.00000
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 15 Rel. abun. bar plot - regular~
Need to make some edits to ‘culture_name’ column for plotting. Specifically, adding wolbachia strain names
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 161 taxa and 458 samples ]:
## sample_data() Sample Data: [ 458 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 161 taxa by 11 taxonomic ranks ]:
## taxa are columnstax.clean <- data.frame(ps.clean@tax_table)
#write.csv(tax.clean,file="tax.clean.toedit.csv")
##read in with edits
tax.edits <- read.csv("tax.clean.edits.csv",row.names=1)
ps.clean.copy <- ps.clean
ps.clean.copy@tax_table <- tax_table(as.matrix(tax.edits))
ps.mq1 <- subset_samples(ps.clean.copy,type=="A.albopictus")
ps.mq <- prune_taxa(taxa_sums(ps.mq1) > 0, ps.mq1)
ps.mq## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 89 taxa and 195 samples ]:
## sample_data() Sample Data: [ 195 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 89 taxa by 11 taxonomic ranks ]:
## taxa are columnsps.mw1 <- subset_samples(ps.clean.copy,type=="Microbial Water")
ps.mw <- prune_taxa(taxa_sums(ps.mw1) > 0, ps.mw1)
ps.mw## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 49 taxa and 211 samples ]:
## sample_data() Sample Data: [ 211 samples by 24 sample variables ]:
## tax_table() Taxonomy Table: [ 49 taxa by 11 taxonomic ranks ]:
## taxa are columns5.1 Prep - mosquitoes
##Having OTU & size columns in front of Kingdom messes up phyloseq stuff
tax.mq <- data.frame(ps.mq@tax_table)
tax.mq.cut <- tax.mq[,3:10]
tax.mq.cut$OTU <- row.names(tax.mq.cut)
##custom grouping for microshades below
tax.mq.cut$gen_culture <- paste0(tax.mq.cut$Genus,"_",tax.mq.cut$Pool_name)
##making a copy of ps object so I don't overwrite original
ps.mq.order <- ps.mq
ps.mq.order@tax_table <- tax_table(as.matrix(tax.mq.cut))
##relative abundance
ps.mq.rel <- transform_sample_counts(ps.mq.order, function(x) x / sum(x))
samdf.mq <- data.frame(ps.mq.rel@sam_data)
samdf.mq$glom <- paste0(samdf.mq$sex,"_",samdf.mq$infusion,"_",samdf.mq$temperature)
samdf.mq$glom## [1] "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C" "M_OL_C" "M_OL_C" "M_OL_C" "F_OL_C"
## [9] "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C"
## [17] "M_OL_C" "M_OL_H" "M_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "M_OL_H" "M_OL_H"
## [25] "F_OL_H" "F_OL_H" "M_OL_H" "F_OL_H" "M_OL_H" "M_OL_C" "F_OL_C" "M_OL_C"
## [33] "M_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C"
## [41] "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_H" "M_OL_H" "F_OL_H" "F_OL_H"
## [49] "F_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "F_OL_H" "F_OL_H" "M_OL_H"
## [57] "M_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "F_OL_H" "M_OL_H"
## [65] "M_OL_H" "F_OL_H" "F_OL_H" "M_OL_C" "F_OL_C" "M_OL_C" "F_OL_C" "F_OL_C"
## [73] "F_OL_C" "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "M_OL_C" "M_OL_C" "M_OL_C"
## [81] "M_OL_C" "M_OL_C" "F_OL_C" "F_OL_C" "F_OL_C" "F_OL_C" "F_OL_H" "M_OL_H"
## [89] "M_OL_H" "F_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "M_OL_H"
## [97] "F_OL_H" "F_OL_H" "M_OL_H" "M_OL_H" "F_OL_H" "M_SG_C" "M_SG_C" "M_SG_H"
## [105] "F_SG_C" "M_SG_C" "F_SG_C" "F_SG_C" "F_SG_C" "F_SG_C" "F_SG_C" "F_SG_C"
## [113] "M_SG_C" "M_SG_C" "M_SG_C" "M_SG_C" "F_SG_C" "F_SG_C" "M_SG_C" "F_SG_C"
## [121] "F_SG_H" "M_SG_H" "M_SG_H" "M_SG_H" "F_SG_H" "F_SG_H" "M_SG_H" "F_SG_H"
## [129] "F_SG_H" "F_SG_H" "M_SG_H" "M_SG_H" "F_SG_H" "F_SG_H" "F_SG_H" "M_SG_H"
## [137] "F_SG_H" "F_SG_H" "M_SG_H" "M_SG_C" "F_SG_H" "F_SG_H" "F_SG_H" "M_SG_H"
## [145] "M_SG_H" "F_SG_H" "F_SW_C" "M_SW_C" "F_SW_C" "M_SW_C" "M_SW_C" "M_SW_C"
## [153] "F_SW_C" "M_SW_C" "M_SW_C" "F_SW_H" "F_SW_H" "F_SW_H" "F_SW_C" "F_SW_C"
## [161] "M_SW_C" "M_SW_C" "F_SW_C" "M_SW_C" "M_SW_C" "F_SW_C" "F_SW_C" "F_SW_C"
## [169] "M_SW_C" "M_SW_C" "M_SW_C" "M_SW_C" "F_SW_C" "M_SW_C" "M_SW_C" "M_SW_C"
## [177] "F_SW_H" "M_SW_H" "F_SW_H" "M_SW_C" "M_SW_C" "M_SW_C" "M_SW_C" "F_SW_C"
## [185] "M_SW_C" "M_SW_C" "F_SW_C" "F_SW_C" "M_SW_C" "F_SW_H" "F_SW_H" "F_SW_H"
## [193] "F_SW_H" "M_SW_H" "F_SW_H"ps.mq.rel@sam_data <- sample_data(samdf.mq)
ps.mq.rel.glom <- merge_samples2(ps.mq.rel, "glom")
ps.mq.rel.glom.rel <- transform_sample_counts(ps.mq.rel.glom, function(x) x / sum(x))
plot_bar(ps.mq.rel.glom.rel,fill="Family")## Warning in psmelt(physeq, as = "data.frame"): The rank names:
## OTU
## have been renamed to:
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.
##Who are the top families for microshades
# Merges ASVs that have the same taxonomy rank (Genus)
ps.mq.rel.fam <- tax_glom(ps.mq.rel, taxrank = "Family") #28 families, 39 gen.
# Calculate taxa sum
top5 = head(sort(colSums(otu_table(ps.mq.rel.fam)), decreasing = TRUE), 10)
# Combine count and taxonomyTable
top5 = cbind(as.data.frame(tax_table(ps.mq.rel.fam)[names(top5),]), Count = top5)
top5## Class Order Family Genus Species
## Otu0001 Alphaproteobacteria Rickettsiales Anaplasmataceae <NA> <NA>
## Otu0003 Gammaproteobacteria Xanthomonadales Xanthomonadaceae <NA> <NA>
## Otu0005 Bacteroidia Flavobacteriales Weeksellaceae <NA> <NA>
## Otu0002 Gammaproteobacteria Enterobacterales Enterobacteriaceae <NA> <NA>
## Otu0022 Bacilli Lactobacillales Carnobacteriaceae <NA> <NA>
## Otu0010 Gammaproteobacteria Pseudomonadales Pseudomonadaceae <NA> <NA>
## Otu0031 Alphaproteobacteria Rhizobiales Xanthobacteraceae <NA> <NA>
## Otu0121 Gammaproteobacteria Enterobacterales Orbaceae <NA> <NA>
## Otu0007 Alphaproteobacteria Acetobacterales Acetobacteraceae <NA> <NA>
## Otu0137 Gammaproteobacteria Burkholderiales Burkholderiaceae <NA> <NA>
## Sequence OTU Pool_name gen_culture Count
## Otu0001 <NA> <NA> <NA> <NA> 170.19659314
## Otu0003 <NA> <NA> <NA> <NA> 11.05234660
## Otu0005 <NA> <NA> <NA> <NA> 5.79876664
## Otu0002 <NA> <NA> <NA> <NA> 3.93453195
## Otu0022 <NA> <NA> <NA> <NA> 1.85376844
## Otu0010 <NA> <NA> <NA> <NA> 1.83582760
## Otu0031 <NA> <NA> <NA> <NA> 0.14403788
## Otu0121 <NA> <NA> <NA> <NA> 0.03718247
## Otu0007 <NA> <NA> <NA> <NA> 0.03185857
## Otu0137 <NA> <NA> <NA> <NA> 0.020614695.2 Microshades - mosquitoes
## Warning in speedyseq::psmelt(.): The rank names:
## OTU
## have been renamed to:
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.# Create a color object for the specified data
col.mdf.mq <- create_color_dfs(mdf.mq, top_orientation = FALSE,group_level="Family",subgroup_level="gen_culture",selected_groups=c("Anaplasmataceae","Xanthomonadaceae","Weeksellaceae","Enterobacteriaceae"),cvd=TRUE)
#Extract
col.mdf.mq.m <- col.mdf.mq$mdf
col.mdf.mq.c <- col.mdf.mq$cdf
##default plot
plot_microshades(col.mdf.mq.m, col.mdf.mq.c)+
#facet_wrap(scales="free")+
theme_cowplot()+
facet_wrap(~infusion,scales="free")+
scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))+
xlab("")+
ggtitle("Mosquitoes")+
guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))
#
# ##sample data got removed during conglomerating
# col.mdf.mq.m$sex <- substr(col.mdf.mq.m$Sample,1,1)
# col.mdf.mq.m$infusion <- substr(col.mdf.mq.m$Sample,3,4)
# col.mdf.mq.m$temperature <- substr(col.mdf.mq.m$Sample,6,6)
# col.mdf.mq.m$infusion <- sub("SW","PW",col.mdf.mq.m$infusion)
# col.mdf.mq.m$infusion <- factor(col.mdf.mq.m$infusion,levels=c("OL","SG","PW"))
col.mdf.mq.c$group## [1] "Enterobacteriaceae-Kosakonia_KOSA1"
## [2] "Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2"
## [3] "Enterobacteriaceae-Klebsiella_KLEB1"
## [4] "Enterobacteriaceae-Family_Enterobacteriaceae_"
## [5] "Enterobacteriaceae-Other"
## [6] "Weeksellaceae-Chryseobacterium_CHRY1"
## [7] "Weeksellaceae-Chryseobacterium_"
## [8] "Xanthomonadaceae-Stenotrophomonas_STEN1"
## [9] "Xanthomonadaceae-Stenotrophomonas_STEN2"
## [10] "Xanthomonadaceae-Stenotrophomonas_"
## [11] "Xanthomonadaceae-Thermomonas_XANT1"
## [12] "Anaplasmataceae-Wolbachia_wAlbB"
## [13] "Anaplasmataceae-Wolbachia_wAlbA"
## [14] "Anaplasmataceae-Wolbachia_"
## [15] "Other-Carnobacterium_CARN1"
## [16] "Other-Pseudomonas_PSEU2"
## [17] "Other-Bradyrhizobium_"
## [18] "Other-Asaia_ASAI1/2"
## [19] "Other-Other"##can't see these ones:
bye <- c("Enterobacteriaceae-Family_Enterobacteriaceae_", "Enterobacteriaceae-Other", "Weeksellaceae-Chryseobacterium_", "Xanthomonadaceae-Stenotrophomonas_", "Xanthomonadaceae-Thermomonas_XANT1", "Anaplasmataceae-Wolbachia_", "Other-Bradyrhizobium_", "Other-Asaia_ASAI1/2", "Other-Other")
col.mdf.mq.m.less <- col.mdf.mq.m[!col.mdf.mq.m$group %in% bye,]
##fix the lengths of these ones:
# col.mdf.mq.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_ENTE1/2","Enterobacteriaceae-unclassified_ENTE1/2",col.mdf.mq.m.less$group)
##desired order
col.mdf.mq.m.less$group <- factor(col.mdf.mq.m.less$group,levels=c("Other-Carnobacterium_CARN1","Other-Pseudomonas_PSEU2","Anaplasmataceae-Wolbachia_wAlbB","Anaplasmataceae-Wolbachia_wAlbA","Xanthomonadaceae-Stenotrophomonas_STEN1","Xanthomonadaceae-Stenotrophomonas_STEN2","Weeksellaceae-Chryseobacterium_CHRY1","Enterobacteriaceae-Kosakonia_KOSA1","Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2","Enterobacteriaceae-Klebsiella_KLEB1"))
##color order
col.mdf.mq.c$hex## [1] "#7D3560" "#A1527F" "#CC79A7" "#E794C1" "#EFB6D6" "#148F77" "#009E73"
## [8] "#098BD9" "#56B4E9" "#7DCCFF" "#BCE1FF" "#9D654C" "#C17754" "#F09163"
## [15] "#616161" "#8B8B8B" "#B7B7B7" "#D6D6D6" "#F5F5F5"hex.colors <- c("#616161","#8B8B8B","#9D654C","#C17754","#098BD9","#56B4E9","#148F77","#7D3560","#A1527F","#CC79A7")
gg.bars.mosq <- ggplot(data=col.mdf.mq.m.less,aes(x=Sample,y=Abundance,fill=group,color=group))+
geom_bar(stat="identity", position="stack")+
facet_wrap(~infusion,scales="free")+
scale_fill_manual(values=hex.colors)+
scale_color_manual(values=hex.colors)+
ggtitle("Mosquitoes")+
guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))+
theme_cowplot()+
scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
xlab("")+
theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))
gg.bars.mosq
5.3 Prep - water
##Having OTU & size columns in front of Kingdom messes up phyloseq stuff
tax.mw <- data.frame(ps.mw@tax_table)
tax.mw.cut <- tax.mw[,3:10]
tax.mw.cut$OTU <- row.names(tax.mw.cut)
##custom grouping for microshades below
tax.mw.cut$gen_culture <- paste0(tax.mw.cut$Genus,"_",tax.mw.cut$Pool_name)
##making a copy of ps object so I don't overwrite original
ps.mw.order <- ps.mw
ps.mw.order@tax_table <- tax_table(as.matrix(tax.mw.cut))
##relative abundance
ps.mw.rel <- transform_sample_counts(ps.mw.order, function(x) x / sum(x))
##new sample type so I can group them
samdf.mw <- data.frame(ps.mw.rel@sam_data)
samdf.mw$glom <- paste0(samdf.mw$day,"_",samdf.mw$infusion,"_",samdf.mw$temperature)
samdf.mw$glom## [1] "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_H" "12_OL_H" "12_OL_H"
## [8] "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_H" "12_OL_H" "12_OL_H"
## [15] "12_OL_H" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_C" "12_OL_H" "12_OL_H"
## [22] "12_OL_H" "12_OL_H" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_H"
## [29] "12_SG_H" "12_SG_H" "12_SG_H" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_C"
## [36] "12_SG_H" "12_SG_H" "12_SG_H" "12_SG_C" "12_SG_C" "12_SG_C" "12_SG_C"
## [43] "12_SG_H" "12_SG_H" "12_SG_H" "12_SG_H" "12_SW_C" "12_SW_C" "12_SW_C"
## [50] "12_SW_C" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_C" "12_SW_C"
## [57] "12_SW_C" "12_SW_C" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_C"
## [64] "12_SW_C" "12_SW_C" "12_SW_C" "12_SW_H" "12_SW_H" "12_SW_H" "12_SW_H"
## [71] "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_H" "20_OL_H" "20_OL_H"
## [78] "20_OL_H" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_H" "20_OL_H"
## [85] "20_OL_H" "20_OL_H" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_C" "20_OL_H"
## [92] "20_OL_H" "20_OL_H" "20_OL_H" "20_SG_C" "20_SG_C" "20_SG_C" "20_SG_C"
## [99] "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_C" "20_SG_C" "20_SG_C"
## [106] "20_SG_C" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_C" "20_SG_C"
## [113] "20_SG_C" "20_SG_C" "20_SG_H" "20_SG_H" "20_SG_H" "20_SG_H" "20_SW_C"
## [120] "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_H" "20_SW_H" "20_SW_H" "20_SW_C"
## [127] "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_H" "20_SW_H" "20_SW_H" "20_SW_H"
## [134] "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_C" "20_SW_H" "20_SW_H" "20_SW_H"
## [141] "20_SW_H" "4_OL_C" "4_OL_C" "4_OL_C" "4_OL_C" "4_OL_H" "4_OL_H"
## [148] "4_OL_H" "4_OL_H" "4_OL_C" "4_OL_C" "4_OL_C" "4_OL_C" "4_OL_H"
## [155] "4_OL_H" "4_OL_H" "4_OL_H" "4_OL_C" "4_OL_C" "4_OL_C" "4_OL_C"
## [162] "4_OL_H" "4_OL_H" "4_OL_H" "4_OL_H" "4_SG_C" "4_SG_C" "4_SG_C"
## [169] "4_SG_C" "4_SG_H" "4_SG_H" "4_SG_H" "4_SG_C" "4_SG_C" "4_SG_C"
## [176] "4_SG_H" "4_SG_H" "4_SG_H" "4_SG_H" "4_SG_C" "4_SG_C" "4_SG_C"
## [183] "4_SG_C" "4_SG_H" "4_SG_H" "4_SG_H" "4_SG_H" "4_SW_C" "4_SW_C"
## [190] "4_SW_C" "4_SW_C" "4_SW_H" "4_SW_H" "4_SW_H" "4_SW_H" "4_SW_C"
## [197] "4_SW_C" "4_SW_C" "4_SW_C" "4_SW_H" "4_SW_H" "4_SW_H" "4_SW_H"
## [204] "4_SW_C" "4_SW_C" "4_SW_C" "4_SW_C" "4_SW_H" "4_SW_H" "4_SW_H"
## [211] "4_SW_H"ps.mw.rel@sam_data <- sample_data(samdf.mw)
ps.mw.rel.glom <- merge_samples2(ps.mw.rel, "glom")
ps.mw.rel.glom.rel <- transform_sample_counts(ps.mw.rel.glom, function(x) x / sum(x))
plot_bar(ps.mw.rel.glom.rel,fill="Family")## Warning in psmelt(physeq, as = "data.frame"): The rank names:
## OTU
## have been renamed to:
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.
##Who are the top families for microshades
# Merges ASVs that have the same taxonomy rank
ps.mw.rel.fam <- tax_glom(ps.mw.rel, taxrank = "Family") #28 families, 39 gen.
# Calculate taxa sum
top5 = head(sort(colSums(otu_table(ps.mw.rel.fam)), decreasing = TRUE), 10)
# Combine count and taxonomyTable
top5 = cbind(as.data.frame(tax_table(ps.mw.rel.fam)[names(top5),]), Count = top5)
top5## Class Order Family Genus
## Otu0002 Gammaproteobacteria Enterobacterales Enterobacteriaceae <NA>
## Otu0003 Gammaproteobacteria Xanthomonadales Xanthomonadaceae <NA>
## Otu0007 Alphaproteobacteria Acetobacterales Acetobacteraceae <NA>
## Otu0005 Bacteroidia Flavobacteriales Weeksellaceae <NA>
## Otu0010 Gammaproteobacteria Pseudomonadales Pseudomonadaceae <NA>
## Otu0014 Gammaproteobacteria Enterobacterales Erwiniaceae <NA>
## Otu0031 Alphaproteobacteria Rhizobiales Xanthobacteraceae <NA>
## Otu0001 Alphaproteobacteria Rickettsiales Anaplasmataceae <NA>
## Otu0070 Gammaproteobacteria Enterobacterales Order_Enterobacterales <NA>
## Otu0162 Gammaproteobacteria Enterobacterales Morganellaceae <NA>
## Species Sequence OTU Pool_name gen_culture Count
## Otu0002 <NA> <NA> <NA> <NA> <NA> 149.4822380
## Otu0003 <NA> <NA> <NA> <NA> <NA> 32.0710953
## Otu0007 <NA> <NA> <NA> <NA> <NA> 13.1670955
## Otu0005 <NA> <NA> <NA> <NA> <NA> 6.4043253
## Otu0010 <NA> <NA> <NA> <NA> <NA> 4.3610484
## Otu0014 <NA> <NA> <NA> <NA> <NA> 3.3758270
## Otu0031 <NA> <NA> <NA> <NA> <NA> 1.0753763
## Otu0001 <NA> <NA> <NA> <NA> <NA> 0.7593486
## Otu0070 <NA> <NA> <NA> <NA> <NA> 0.2676740
## Otu0162 <NA> <NA> <NA> <NA> <NA> 0.02557955.4 Microshades - water
## Warning in speedyseq::psmelt(.): The rank names:
## OTU
## have been renamed to:
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.# Create a color object for the specified data
col.mdf.mw <- create_color_dfs(mdf.mw, top_orientation = FALSE,group_level="Family",subgroup_level="gen_culture",selected_groups=c("Enterobacteriaceae","Xanthomonadaceae","Acetobacteraceae","Weeksellaceae"),cvd=TRUE)
#Extract
col.mdf.mw.m <- col.mdf.mw$mdf
col.mdf.mw.c <- col.mdf.mw$cdf
##default plot
plot_microshades(col.mdf.mw.m, col.mdf.mw.c)+
#facet_wrap(scales="free")+
theme_cowplot()+
#scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))+
xlab("")+
guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))
col.mdf.mw.c.new <- color_reassign(col.mdf.mw.c,
group_assignment = c("Enterobacteriaceae","Xanthomonadaceae","Acetobacteraceae","Weeksellaceae"),
color_assignment = c("micro_cvd_purple","micro_cvd_blue","micro_cvd_green","micro_cvd_turquoise"),group_level="Family")
##check new colors
plot_microshades(col.mdf.mw.m, col.mdf.mw.c.new)+
#facet_wrap(scales="free")+
theme_cowplot()+
#scale_x_discrete(labels=c("F_cool","F_warm","M_cool","M_warm"))+
theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))+
xlab("")+
guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))
# ##sample data got removed during conglomerating
# col.mdf.mw.m2 <- cbind(col.mdf.mw.m,data.frame(str_split_fixed(col.mdf.mw.m$Sample,"_",3)))
# col.mdf.mw.m2$X1==col.mdf.mw.m2$day #day lines up
# #X2 = infusion, X3=temp
col.mdf.mw.m$infusion <- sub("SW","PW",col.mdf.mw.m$infusion)
col.mdf.mw.m$infusion <- factor(col.mdf.mw.m$infusion,levels=c("OL","SG","PW"))
col.mdf.mw.c.new$group## [1] "Weeksellaceae-Chryseobacterium_CHRY1"
## [2] "Acetobacteraceae-Asaia_ASAI1/2"
## [3] "Acetobacteraceae-Asaia_ASAI3"
## [4] "Acetobacteraceae-Asaia_"
## [5] "Xanthomonadaceae-Stenotrophomonas_STEN1"
## [6] "Xanthomonadaceae-Stenotrophomonas_STEN2"
## [7] "Xanthomonadaceae-Stenotrophomonas_"
## [8] "Xanthomonadaceae-Thermomonas_XANT1"
## [9] "Enterobacteriaceae-Kosakonia_KOSA1"
## [10] "Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2"
## [11] "Enterobacteriaceae-Kosakonia_KOSA2"
## [12] "Enterobacteriaceae-Klebsiella_KLEB1"
## [13] "Enterobacteriaceae-Other"
## [14] "Other-Pantoea_PANT1/2/3"
## [15] "Other-Pseudomonas_PSEU2"
## [16] "Other-Bradyrhizobium_"
## [17] "Other-Pseudomonas_"
## [18] "Other-Other"##can't see these ones:
bye.mw <- c("Acetobacteraceae-Asaia_","Xanthomonadaceae-Stenotrophomonas_","Xanthomonadaceae-Thermomonas_XANT1","Other-Pseudomonas_","Other-Other")
col.mdf.mw.m.less <- col.mdf.mw.m[!col.mdf.mw.m$group %in% bye.mw,]
# ##fix the lengths of these ones:
# ##Entero group
# col.mdf.mw.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_ENTE1/2","Enterobacteriaceae-unclassified_ENTE1/2",col.mdf.mw.m.less$group)
# col.mdf.mw.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_KLEB1","Enterobacteriaceae-Klebsiella_KLEB1",col.mdf.mw.m.less$group)
# col.mdf.mw.m.less$group <- sub("Enterobacteriaceae-Enterobacteriaceae_unclassified_KOSA2","Enterobacteriaceae-Kosakonia_KOSA2",col.mdf.mw.m.less$group)
##Other group
col.mdf.mw.m.less$group <- sub("Other-Bradyrhizobium_","Other-Bradyrhizobium",col.mdf.mw.m.less$group)
# col.mdf.mw.m.less$group <- sub("Other-Gammaproteobacteria_unclassified_","Other-Gammaproteobacteria_unclassified",col.mdf.mw.m.less$group)
##desired order
col.mdf.mw.m.less$group <- factor(col.mdf.mw.m.less$group,levels=c("Other-Pantoea_PANT1/2/3","Other-Pseudomonas_PSEU2","Other-Bradyrhizobium","Acetobacteraceae-Asaia_ASAI1/2","Acetobacteraceae-Asaia_ASAI3","Xanthomonadaceae-Stenotrophomonas_STEN1","Xanthomonadaceae-Stenotrophomonas_STEN2","Weeksellaceae-Chryseobacterium_CHRY1","Enterobacteriaceae-Kosakonia_KOSA1","Enterobacteriaceae-Family_Enterobacteriaceae_ENTE1/2","Enterobacteriaceae-Klebsiella_KLEB1","Enterobacteriaceae-Kosakonia_KOSA2","Enterobacteriaceae-Other"))
##switching some colors to match mosquitoes
hex.colors.mw <- c("#616161","#8B8B8B","#B7B7B7","#4E7705","#6D9F06","#098BD9","#56B4E9","#148F77","#7D3560","#A1527F","#CC79A7","#E794C1","#EFB6D6")
##Sample order
col.mdf.mw.m.less$Sample <- factor(col.mdf.mw.m.less$Sample, levels=c("4_OL_C","4_OL_H","12_OL_C","12_OL_H","20_OL_C","20_OL_H","4_SG_C","4_SG_H","12_SG_C","12_SG_H","20_SG_C","20_SG_H","4_SW_C","4_SW_H","12_SW_C","12_SW_H","20_SW_C","20_SW_H"))
gg.bars.water <- ggplot(data=col.mdf.mw.m.less,aes(x=Sample,y=Abundance,fill=group,color=group))+
geom_bar(stat="identity", position="stack")+
facet_wrap(~infusion,scales="free")+
scale_fill_manual(values=hex.colors.mw)+
scale_color_manual(values=hex.colors.mw)+
ggtitle("Mesocosm water")+
guides(color=guide_legend("Family-Genus-Culture"),fill=guide_legend("Family-Genus-Culture"))+
theme_cowplot()+
scale_x_discrete(labels=c("D4_cool","D4_warm","D12_cool","D12_warm","D20_cool","D20_warm"))+
xlab("")+
theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))
gg.bars.water
6 Rel. abun. box plots - mosquitoes
ps.trim.mq.rel <- transform_sample_counts(ps.trim.mq, function(x) x / sum(x))
melt.trim.mq.rel <- psmelt(ps.trim.mq.rel)## Warning in psmelt(ps.trim.mq.rel): The rank names:
## OTU
## have been renamed to:
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.6.1 Wolb in sex
melt.trim.mq.rel.wolb <- subset(melt.trim.mq.rel,Genus=="Wolbachia")
ggplot(data = melt.trim.mq.rel.wolb, aes(x = sex, y = Abundance))+
geom_boxplot(outlier.shape = NA)+
geom_jitter(aes(color = OTU), height = 0, width = .2) +
facet_wrap(~OTU,scales="free")
melt.trim.mq.rel.1.4 <- subset(melt.trim.mq.rel, OTU %in% c("Otu0001", "Otu0004"))
melt.trim.mq.rel.1.4$Longer_name <- melt.trim.mq.rel.1.4$OTU
melt.trim.mq.rel.1.4$Longer_name <- sub("Otu0001","Otu0001_wAlbB",melt.trim.mq.rel.1.4$Longer_name)
melt.trim.mq.rel.1.4$Longer_name <- sub("Otu0004","Otu0004_wAlbA",melt.trim.mq.rel.1.4$Longer_name)
gg.sex.wolb <- ggplot(data = melt.trim.mq.rel.1.4, aes(x = sex, y = Abundance))+
geom_jitter(width=0.2,alpha=0.5)+
geom_boxplot(outlier.shape = NA,alpha=0.5)+
facet_wrap(~Longer_name,scales="free")+
xlab("Mosquito sex")+
ylab("Relative abundance/sample")+
theme_bw()+
scale_x_discrete(labels=c("Female","Male"))+
guides(color="none")
gg.sex.wolb
6.2 Chry & pseu in infusion
melt.trim.mq.rel.5.10 <- subset(melt.trim.mq.rel, OTU %in% c("Otu0005", "Otu0010"))
gg.inf <- ggplot(data = melt.trim.mq.rel.5.10, aes(x = infusion, y = Abundance))+
geom_jitter(width=0.2,alpha=0.5)+
geom_boxplot(outlier.shape = NA,alpha=0.5)+
facet_wrap(~Longer_name,scales="free")+
ylab("Relative abundance/sample")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
theme_bw()+
guides(color="none")
gg.inf
6.3 Kleb in temp
melt.trim.mq.rel.kleb <- subset(melt.trim.mq.rel, OTU %in% c("Otu0009","Otu0022"))
gg.tem <- ggplot(data = melt.trim.mq.rel.kleb, aes(x = temperature, y = Abundance))+
geom_jitter(width=0.2,alpha=0.5)+
geom_boxplot(outlier.shape = NA,alpha=0.5)+
facet_wrap(~Longer_name,scales="free")+
ylab("Relative abundance/sample")+
theme_bw()+
xlab("Temperature")+
scale_x_discrete(labels=c("Cool","Warm"))+
ylim(0,0.017)+
guides(color="none")
gg.tem ## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 182 rows containing missing values or values outside the scale range
## (`geom_point()`).
7 Relative abundance box plots for mesocosm water
ps.trim.mw.rel <- transform_sample_counts(ps.trim.mw, function(x) x / sum(x))
melt.trim.mw.rel <- psmelt(ps.trim.mw.rel)## Warning in psmelt(ps.trim.mw.rel): The rank names:
## OTU
## have been renamed to:
## taxa_OTU
## to avoid conflicts with special phyloseq plot attribute names.7.1 Infusion differences
melt.trim.mw.inf <- subset(melt.trim.mw.rel, OTU %in% c("Otu0003","Otu0005","Otu0007","Otu0024"))
ggplot(data = melt.trim.mw.inf, aes(x = infusion, y = Abundance))+
geom_boxplot(outlier.shape = NA)+
geom_jitter(aes(color = Longer_name), height = 0, width = .2) +
facet_wrap(~Longer_name,scales="free")
gg.mw.inf <- ggplot(data = melt.trim.mw.inf, aes(x = infusion, y = Abundance))+
geom_jitter(width=0.2,alpha=0.5)+
geom_boxplot(outlier.shape = NA,alpha=0.5)+
facet_wrap(~Longer_name,scales="free")+
xlab("Infusion")+
ylab("Relative abundance/sample")+
theme_bw()+
scale_x_discrete(labels=c("OL","SG","PW"))+
guides(color="none")
gg.mw.inf
8 Session info
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.2.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Pacific/Honolulu
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] speedyseq_0.5.3.9021 stringr_1.5.1 microshades_1.13
## [4] ggpubr_0.6.0 microViz_0.12.4 dplyr_1.1.4
## [7] cowplot_1.1.3 ggplot2_3.5.1 phyloseq_1.48.0
## [10] vegan_2.6-6.1 lattice_0.22-6 permute_0.9-7
##
## loaded via a namespace (and not attached):
## [1] ade4_1.7-22 tidyselect_1.2.1 farver_2.1.2
## [4] Biostrings_2.72.1 fastmap_1.2.0 digest_0.6.37
## [7] lifecycle_1.0.4 cluster_2.1.6 survival_3.7-0
## [10] magrittr_2.0.3 compiler_4.4.1 rlang_1.1.4
## [13] sass_0.4.9 tools_4.4.1 igraph_2.0.3
## [16] utf8_1.2.4 yaml_2.3.10 data.table_1.15.4
## [19] knitr_1.48 ggsignif_0.6.4 labeling_0.4.3
## [22] plyr_1.8.9 abind_1.4-5 Rtsne_0.17
## [25] withr_3.0.1 purrr_1.0.2 BiocGenerics_0.50.0
## [28] grid_4.4.1 stats4_4.4.1 fansi_1.0.6
## [31] multtest_2.60.0 biomformat_1.32.0 colorspace_2.1-1
## [34] Rhdf5lib_1.26.0 scales_1.3.0 iterators_1.0.14
## [37] MASS_7.3-61 cli_3.6.3 rmarkdown_2.28
## [40] crayon_1.5.3 generics_0.1.3 microbiome_1.26.0
## [43] httr_1.4.7 reshape2_1.4.4 ape_5.8
## [46] cachem_1.1.0 rhdf5_2.48.0 zlibbioc_1.50.0
## [49] splines_4.4.1 parallel_4.4.1 XVector_0.44.0
## [52] vctrs_0.6.5 Matrix_1.7-0 carData_3.0-5
## [55] jsonlite_1.8.8 car_3.1-2 IRanges_2.38.1
## [58] S4Vectors_0.42.1 rstatix_0.7.2 foreach_1.5.2
## [61] tidyr_1.3.1 jquerylib_0.1.4 glue_1.7.0
## [64] codetools_0.2-20 stringi_1.8.4 gtable_0.3.5
## [67] GenomeInfoDb_1.40.1 UCSC.utils_1.0.0 munsell_0.5.1
## [70] tibble_3.2.1 pillar_1.9.0 htmltools_0.5.8.1
## [73] rhdf5filters_1.16.0 GenomeInfoDbData_1.2.12 R6_2.5.1
## [76] evaluate_0.24.0 Biobase_2.64.0 highr_0.11
## [79] backports_1.5.0 broom_1.0.6 bslib_0.8.0
## [82] Rcpp_1.0.13 gridExtra_2.3 nlme_3.1-166
## [85] mgcv_1.9-1 xfun_0.47 forcats_1.0.0
## [88] pkgconfig_2.0.3